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Abstract 

A new pseudospectral calculation of collapsing Euler vortices Hou &: Li (2006) has 
called into question the long-term conclusions of singular behavior described earlier in 
Kerr (1993, 2005). This review is designed to: improve the discussion of the detailed 
analysis of one test initial condition designed find sources of errors, to compare that 
with calculations showing no evidence of a singularity, and to document two sets of 
discussions. Those prior to 1993 between competing teams. And recent discussions of 
what is needed to reach more convincing conclusions. 
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' — ^1 Introduction 

J> Whether or not there are singularities of the incompressible three-dimensional Euler equa- 
^ tions is of fundamental mathematical and physical importance. This review will not discuss 
^ Either of these issues or related questions for how this behavior could be applied to our under- 
^ standing of partial differential equations or turbulence theory. Instead, the focus will be the 
\Q history, circa 1990, of numerical calculations that were designed to determine if there could 
O be a singularity of the three-dimensional incompressible Euler equations and how, through 

^ several workshops, a concensus was reached about to address the issue numerically. Analytic 

• ^ results will not be discussed other than how they pertain to these calculations. One test 
J>^talculation from Kerr (1992) will be used to suggest that the differences between the struc- 

^ tures in the new calculation by Hou & Li (2006) and those presented by Kerr (1993) could 
pbe associated with small-scale noise. How this could disrupt a singular structure forming is 

• discussed. 

■ The paper is organized as follows. After the early history of Euler calculations is given, 
a hiethodologies agreed upon by concensus at several workshops are presented, including how 
to compare calculations. Next is a discussion of the uses and dangers of the most robust 
approach, namely pseudospectral methods. Then a section compares Kerr (1993) to compet- 
itive calculations from the early 1990s. Finally there is a comparison between Kerr (1993) 
and Hou & Li (2006). 



2 Early Euler calculations 

The first serious attempt to study Euler numerically was the Taylor-Green calculations of 
Brachet et al. (1983). Their conclusion was that vortex sheets suppressed any trends to- 
wards singularities that had previously been suggested either by spectral closures or series 
expansions. This conclusion was qualitatively similar to the conclusion from two-dimensional 



ideal, incompressible MHD calculations (Sulem et al. (1985)) where it was found that cur- 
rent sheets suppressed trends towards singularities. 

The next set of numerical investigations used vortex filaments. While these are a poor 
approximation of the Euler equations at the smallest scales, they did provide a global pic- 
ture that an anti-parallel initial condition might predispose a calculation towards singular 
behavior in a manner that smooth initial conditions might not (Pumir & Siggia (1987)). 
This led to two preliminary studies at low resolution (Ashurst & Mciron (1987); Pumir 
& Kerr (1987)) whose only significant contribution was that vortex flattening might play 
a role similar to current sheets in 2D MHD. Simultaneously, an orthogonal initial condition 
was simulated that showed more curvature development than the anti-parallel simulations 
(preliminary results by Melander and Zabusky in 1988 eventually appeared in Boratav et al. 

(1992) ). 

All the subsequent calculations assumed an anti-parallel geometry, for which there are 
two symmetry planes. One in y— 2; is between the vortices and was called the 'dividing plane'. 
The other in a; — z is at the position of maximum perturbation and was called the 'symmetry 
plane'. (Terminology due to F. Hussain.) The next step was to find a better initial condition 
within this geometry and the first steps towards adaptivity. Melander & Hussain (1989) 
were able to identify an initial condition where all derivatives in the profile of the vortex 
core went smoothly to zero and used a sinusoidal perturbation of the x — z position as a 
function of y. The calculation did not take advantage of the symmetries, had uniform mesh 
spacing in all direction and was viscous. Kerr & Hussain (1989) tried a truncated Gaussian 
core and resolved the discontinuities at the edge of the initial condition by applying a strong 
Fourier based hyper-viscous filter. Kerr & Hussain (1989) used a combination of higher 
order trigonometric functions to create an initial perturbation that was localized near the 
symmetry plane. Symmetries were used to reduce the computational requirements, the mesh 
spacing depended upon direction and the calculations were viscous. Pumir & Siggia (1990) 
uses a hyperbolic trajectory for the vortex lines, therefore being the only calculation to date 
that is not in a periodic geometry. They used a strictly Gaussian core that does not ever 
go to zero, which is possible if one does not have to worry about overlapping across periodic 
boundaries. Their calculation was completely adaptive and inviscid. 

Of these calculations, only the viscous calculation of Kerr & Hussain (1989) had any con- 
sistency with singular behavior. To test the Melander & Hussain (1989) profile in adaptive 
calculations, Kerr (1993) and Shelley et al. (1993) tried Chebyshev and mapped Fourier 
methods respectively. Kerr (1993) reported consistency with the singular trends of Kerr 
& Hussain (1989). Shelley et al. (1993) concluded that the trends did not favor a singu- 
larity, but their calculations used only the sinusoidal perturbation of Melander & Hussain 
(1989) and were viscous. Why they reached different conclusions was discussed briefly by 
Kerr (1992) and will be discussed further here in subsection 5.3. 

The only significant calculations since then are due to Graucr et al. (1998) and now Hou 
& Li (2006). The Grauer et al. (1998) calculation seems consistent with the trends in Kerr 

(1993) while Hou & Li (2006) is not. Kerr (1993) and Hou & Li (2006) are compared 
in section 6. The Kida vortex promoted by Pelz (2001) as a scenario for a singularity has 
recently been shown to be regular by Grauer (private communication). 

3 Workshop methodologies 

The methodology of Pumir & Siggia (1990) and Kerr (1993) was worked out at two 
workshops on Topological Fluid Dynamics chaired by H. K. Moffatt in 1989 and 1991. These 



were the lUTAM Symposium on Topological Fluid Dynamics in August 1989 in Cambridge, 
England and the Program on Topological Fluid Dynamics at the Institute for Theoretical 
Physics in Santa Barbara in the Fall of 1991, with a symposium on issues in Euler at the 
end of October 1991. 

The following principles were worked out at those programs in discussions between most 
of the principal authors at that time such as U. Frisch, F. Hussain, R.M. Kerr, R. Pelz, A. 
Pumir E.D. Siggia, N. Zabusky, and the other participants in those workshops. 

These are some of the rules: 

a) Run only Euler. Do not try to reach conclusions about Euler using a series of decreas- 
ing viscosity Navier-Stokes calculations. The range of Reynolds numbers available to 
Navier-Stokes calculations at that time is too short to be useful in reaching any definite 
conclusions. 

b) Use refined meshes. Even the easy solution of different mesh sizes in different directions, 
as in Kerr & Hussain (1989), has limitations. There is too much space over which 
nothing is happening. Complementary pseudo-spectral calculations can still be useful 
to confirm the numerical method, but serious compromises have to be made as discussed 
below. 

c) Have a localized initial perturbation. Orthogonal vortices and the hyperbolic trajectory 
of Pumir & Siggia (1990) automatically satisfy this. Kerr & Hussain (1989) and Kerr 
(1993) satisfy this condition using a map of a sinusoidal perturbation. 

d) One needs to demonstrate structures that wouldn't be subject to depletion of nonlin- 
earity. In particular if vortex sheets appear, they must develop strong curvature or 
kinks. 

e) The primary quantity of interest is the norm of vorticity ||cc;||oo so that consistency 
with the analytic result of Beale et al. (1984) can be tested. If singular behavior is 
suspected, then the following must be obeyed: 

t 

\\u>\\oo dt ^ OO (1) 

f) If it has power- law behavior of the form ||a;||oo ~ {T — t)~'^ and 7 < 1, then the 
calculation is incorrect. 

g) The best way to test (1) is to assume 7 = 1. This would be dimensionally consistent. 

h) However, using ||a;||oo by itself can lead to misleading conclusions either way. Double 
exponentials are indistinguishable from power laws and at late times it is impossible to 
avoid some slowing of singular trends, which could lead to the appearance of saturation 
if one assumes that a — d\\u)\\oo /dt. 

i) Therefore one must calculate the time derivative of ||a;||oo directly, that is a = uJieijUJj/\uj 
directly, which should also have time-integral singular growth (Ponce (1985)). In the 
anti-parallel case a is just the velocity derivative on the symmetry plane of the axial 
velocity. This gives an independent measure of singular activity While finding at 
the point of ||a;||oo is most desirable, finding = sup (a) in the symmetry plane was 
found to work better in Kerr (1993). 



j) One needs a measure of collapse to satisfy the condition of Caffarelle et al. (1982). 



Some of these rules were refined by further conversations between R.M. Kerr and A. 
Majda during a two week workshop at the Research Institute in Mathematical Sciences in 
Kyoto, Japan in October 1992. 

Besides the rules just listed, other tests that have been tried and their order of success 
are as follows: 

k) Finding independently that 

with the same singular time as for ||a;||oo and a. The evidence from Kerr (1993) is 
that this test works better than finding either ap or as defined above. 

1) Pumir & Siggia (1990) discusses how the position of ||u;||oo moves. Kerr (1993) 
goes further and tries to show that the positions {x, z)p of all minima and maxima 
of components of the stress tensor in the symmetry plane collapse to the positions of 
||a;||oo ■ That is 

Xp-X{T)-^T-t , Zpr^T-t 
m) Kerr (2005) goes further still by looking at profiles. 

n) sup(|f p) T — t where v is the axial velocity in the direction of vorticity in the 
symmetry plane. This would provide consistency with a suggestion from Constantin et 
al. (1996). The only evidence for this is marginal from Kerr (2005). We need a good 
calculation to test this. 

o) Curvature blowup as k""^ ~ {T — t). Again, this would provide consistency with a 
suggestion from Constantin et al. (1996). So far this is only inferred by related scaling 
properties in Kerr (2005). Again, a good calculation is needed to test this. 



4 Uses and pitfalls for pseudospectral 

For any problem where most of the domain is not involved in the dynamics, as for this 
possibly singular equation, uniform mesh methods waste most of the computational domain. 
Pseudospectral codes fall into this class of methods. Besides being limited in the local reso- 
lution they can provide, pseudospectral codes have the additional annoyance of a continuing 
debate over how to handle the high wavenumbcr cutoff. 

Many strategies have been tried over the years. Initially Orszag and Patterson decided 
upon a spherical truncation which was corrected by adding together two phase-shifted ver- 
sions of the nonlinear terms. Then NASA Ames proposed the Leonard filter for LES cal- 
culations. By the mid-1980s it was generally agreed that none of these tricks gained one 
anything. LES calculations started using only a sharp- wavenumber cutoff with the 2/3rds 
rule and DNS in wall-bounded flows at NASA Ames/CTR did the same. All of my calcula- 
tions after 1985, whether Euler, Navier-Stokes, or convection have followed this pohcy. 

The issue has been reopened by Hou & Li (2006) claiming to have found a filter with 
36th-order accuracy. So let me summarize my understanding of the good and bad points of 
different methods for the Euler problem. 

Options for high wavenumber cutoff of Euler: 



• 2/3rds dealiasing. 

— Good: absolutely no wavenumber contributions for k > {2/3)kmax, which would 
immediately give errors. 

— Good: energy is exactly conserved. 

— Bad: An abrupt cutoff leads to Gibbs phenomena (oscillations) in physical space. 

— Solution: Must set a strict standard when the calculations must be shut off based 
upon anomalous growth of the wavenumbers near the cutoff. This is good in that 
there is a measure of the errors. 

• A smooth filter that ends for k = {2/3)kmax- 

— Good: minimize Gibbs phenomena. 

— Bad: Dissipative. 

— Bad: No measure of errors. 

• A smooth filter that extends beyond k = {2/3)kmax- 

— Good: Maybe minimizes Gibbs phenomena. 

— Bad: Immediately adds aliasing errors. 

— Bad: might add small regions of negative vorticity that can blow up. 

— Bad: No measure of errors. 

— Unsolved: What of the claim in Hou & Li (2006) that this converges? And how 
does one identify what the aliasing errors turn into? 

Despite these problems, it was concluded that in the early stages a pseudo-spectral code is 
much better for testing initial conditions, in particular testing schemes for filtering the initial 
condition. It was in ensuring that the nonlinearity was not depleted even at low resolution 
and early times in such test calculations that led to the initial condition and numerics that 
worked so well in Kerr (1993). 

5 Comparison of structures and ||ci;||oo between older 
calculations 

This section will compare three calculations that meet the simulation criteria (a-c) above: 

• Kerr (1993) 

• An intermediate time from Pumir & Siggia (1990). 

• An inviscid version of the Shelley et al. (1993) calculation. 

Earlier calculations such as Ashurst & Meiron (1987), Pumir & Kerr (1987) and Me- 
lander & Hussain (1989) met none of these criteria on coarse meshes. 

When Kerr (1993) was published, the weight of numerical evidence from simulations with 
similar initial conditions was against a singularity. Therefore, in making claims of singular 
behavior, besides demonstrating points (a-k) above,the following was addressed: 
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Figure 1: ||u;||oo and a contours in the symmetry plane for t = 15 from Kerr (1993) 
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Figure 2: Dependence of l/||a;||oo , 0.1/|ej^j^p| and 0.0003/f2pr upon time from the anti- 
parallel Euler calculation showing convergence to a singular time of about T = 18.7. 




Figure 3: Initial 2;-Chebyshev distribution 
(spectrum) for filtered intial conditions. 



Figure 4: Dependence of ||c(;||oo , \eyy,p\~'^ and 
\^yy,p\ upon time from the filtered initial con- 
dition anti-parallel Euler calculation 



• The computational power of the calculations finding exponential growth was matched. 

• A modification of the singular initial condition was found that was able to reproduce 
the exponential behavior of competitive calculations. 

• A proposal was made for why there could be at least two types of behavior, exponential 
and singular. 

While Pumir & Siggia (1990) and Kerr (1993) meet all three requirements, the published 
work of Shelley et al. (1993) was neither inviscid nor did it have a locahzed perturbation. 
A better description of Shelley et al. (1993) is that it was a shghtly adaptive and higher 
Reynolds number version of Mclandcr & Hussain (1989). Nonetheless, by giving the profile 
they used a localized perturbation and running it inviscidly, the primary conclusion of Shelley 
et al. (1993) could be verified. This was that for the Melander & Hussain (1989) profile, 
singular behavior does not appear. 

The goal of working by steps through these three cases is to suggest how the calculation 
of Hou & Li (2006) might be suppressing singular behavior due to numerical noise in their 
numerical method. The first step follows Kerr (1992) to show how the initial condition 
of Shelley et al. (1993) leads to false regular behavior. Then the graphics for an inviscid 
version of Shelley et al. (1993) is compared to the very late times of Pumir & Siggia (1990). 
In this way the structures associated with regular behavior will be identified. Finally, similar 
structures are identified in Hou & Li (2006). This raises the possibility that the numerical 
method of Hou & Li (2006) also introduces small scale numerical errors. This will only be 
confirmed if their 2/3rds dealiasing calculations are report or are repeated by someone else. 




Figure 5: uj in the symmetry plane, z is not stretched, for i = 0, 6, 12 and 15.5 



Figure 6: ||w||oo and a contours in the symmetry plane for t — 7.2 from Pumir & Siggia 
(1990) 



0.14 
0.12 

o 0.08 

0.06 

0.04 

0.02 
0.00 

6 



G 



1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 








(20cf)-' : 




>, _ 
^ \ 

, 1 . , . . 1 . . , . 1 . . , . 



5,5 7.0 7,5 8,0 8.5 

time 



Figure 7: Growth of ||w||oo , oc from Pumir & Siggia (1990) 



5.1 Comparing structures of Kerr (1993) and Pumir & Siggia 
(1990) 

The purpose of this subsection is to compare the structures that arise during the period 
with the strongest evidence for singular growth of Kerr (1993) with a brief period of growth 
from Pumir & Siggia (1990) that might have consistent behavior. Figure 2 from Kerr 
(2005) shows the growth of ||c(;||oo , iiiax(a) = ||ej,.y|| in the symmetry plane, and the enstrophy 
production flp^- This was the primary evidence that these calculations have singular growth. 
Figure 4 plots the peak vorticity, peak strain along the vorticity in the symmetry plane, and 
the inverse of that strain as functions of time for a 256 x 256 x 192 pseudospectral calculation 
using the initial condition of Kerr (1993). This pseudospectral calculation was not the best 
that was done in preparation for the final calculations of Kerr (1993), but it can be used to 
illustrate the differences between this initial condition and the Melander & Hussain (1989) 
initial condition. The peak of the ctoo = ||eyy||oo strain is plotted instead of the value at the 
position of Halloo because it was better behaved as previously discussed (Kerr (1993)). 

The vorticity and strain contours in Fig. 1 are straight out of Kerr (1993) for t = 15. 
These plots have been expanded by a factor of 5 in the vertical (z) direction. This time was 
chosen because it shows the least numerical noise while being the regime showing singular 
growth of 1 1 I loo • Subsequent to 1993 it was realized that most of the noise in the contours 
in Kerr (1993) comes from a high wavenumber spectral tail, which when filtered out gives 
the smoother figures of Kerr (2005). It is likely that the physical space data for these fields 
in this plane is still stored on the NCAR mass store and might be used later to create better 
figures. 

The two points of comparison to be made are that: 

a) The highest contour of vorticity is placed in the corner of the bent front. 

b) There is a steepening of the positive strain contours in the direction of propagation of 
the front. 

To compare with competitive calculations whose graphic results do not stretch the inter- 
vortical (z) direction, it is useful to show contour plots of vorticity in the symmetry plane 
without stretching in Fig. 5. The evolution with time is as follows: 

• At i = 6 there is only flattening. 

• At t = 12 a tilt in the innermost contour appears. 

• This becomes more of a crescent in the bend at t = 15.5, which is essentially an 
unstretched version of left frame of Fig. 1 with the highest vorticity contour only 
within the bend. 

How does this compare to Pumir & Siggia (1990)? Recall that they use hyperbohc 
trajectories, while all the other calculations use periodic trajectories. Therefore, definitive 
comparisons are out of the question. Nonetheless, Fig. 7 shows that growth of ||a;||oo and 
dip, the strain at the position of ||w||oo , in Pumir & Siggia (1990) is similar to Fig. 2. 
Recall that in Kerr (1993) it was found that while Hej^yHoo ~ 0.1||a;||oo , <^p ~ 0.05||a>||oo ■ 
Therefore there is a short period between t = 7.6 and t = 8.1 where the growth of ||t<^||oo 
and Op in Pumir & Siggia (1990) might be consistent with Kerr (1993). Most of the growth 
of ||u;||oo before t = 7.6 is just to set up the structure, t > 8.1, when saturation appears, is 
discussed in the next subsection. 



Figure 6 shows contours of ||a>||oo and a just before this period sA, t — 7.2. As in Fig. 
1, the highest contour of vorticity is placed in the corner of the bent front and there is a 
steepening of the positive strain contours in the direction of propagation of the front. There 
is a long tail in the vorticity plot, but it seems to be separating from the contour in the 
corner. In these respects the figures are similar to Kerr (1993). 

The next subsection will show that when non-singular growth is observed, these properties 
are not observed. 



5.2 Initial profile and why a localized perturbation with filtering 
is used 

There are two critical components in forming the initial condition previously discussed by 
Kerr (1992) and Kerr (1993). 

• First is using a vortex core profile that, before smoothing, is analytic as originally 
proposed by Melander & Hussain (1989). As they claimed, their profile is far superior 
to the cut-off Gaussian profile used by Kerr & Hussain (1989). However, even here, 
based on spectra (see below), it was decided that this it was not smooth enough and 
a small high-wavenumber filter was necessary. 

• Second, there is the intertwined choices of the initial perturbation in the trajectory of 
the vortex tube and of the dimensions of the periodic domain in the axial direction. 
This is discussed next. 

Following Kerr & Hussain (1989) the trajectory of the center of the initial vortex tubes 
was defined as 

X(s) = Xo + 5x cos(s) 

Z[s) = Zo + 5z cos(s), where 

s{Y) = 1/2 + Ly6yi sin(7ry2/ Ly) and 

y2^Y ^ LySy2 sm{TTY/Ly) 

rather than a simple sinusoidal trajectory. 

The incompressible vorticity vector for this trajectory was proportional to 
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Following a suggestion by Melander (private communication) the profile of vorticity about 
this trajectory was given by 

uj{r) = exp(/(r)). 



where 



fir) 



+ r^(l + r^ 



and a distance r from the center of the vortex core {X, Y, Z) of characteristic radius R is 



r=\{x,y,z)-{X,Y,Z)\/R for r<l 



These formulae have been corrected based upon errors in the text of Kerr (1993) noted 
by Hou & Li (2006). Hon & Li (2006) missed the following additional misprint. In Kerr 
(1993) and Hou & Li (2006) in the formulae for X[s) and Z{s) (or their equivalents) , '/I/^;' 
and '/Lx' respectively should both be ^ / Ly . 

For 5yi = 5y2 = Q the trajectory is sinusoidal, as used by Melander & Hussain (1989). 
Two problems were then identified with the procedure of Melander & Hussain (1989), which 
were replicated by Shelley et al. (1993). 

The first problem identified by Kerr & Hussain (1989) was that if a sinusoidal pertur- 
bation is used then, as the two vortices collapse in the direction between themselves, they 
also expanded in the axial y direction. When this expansion runs into its periodic image, 
growth in II UJ I loo is suppressed. Doubling the domain size in y does not resolve the problem. 
Therefore to localize the perturbation, Kerr & Hussain (1989) and Kerr (1993) chose the 
following values of 5: 5yi = 0.5, 6y2 = 0.4, 6^ = — 1.6and 6^ = and Zo = 1.57 and R = 0.75. 

A more serious problem was that the initial profile of Melander & Hussain (1989), 
without any filtering, created regions of negative vorticity in the symmetry plane. There are 
small negative regions that cannot be seen in the contour plots unless either zero contours 
are plotted or highs and lows are given. This is demonstrated in the first frame of Fig. 8. 

When highs and lows as in Fig 8 are not used in the graphics, one cannot assess the 
possible impact of small negative regions. This would apply to the graphics in Pumir & Siggia 
(1990), Shelley et al. (1993) and now Hou & Li (2006) In Hou & Li (2006) negative 
contours would not appear in the initial condition because their initial condition is filtered, 
but evidence for negative regions at late times is discussed in the next section. 

In preparation for the calculation in Kerr (1993) that used a filtered initial condition, 
a variety of initial conditions without filtering were experimented with. A common feature 
of all the initial z spectra, or equivalently the initial distribution of Chebyshev polynomials, 
was that they went to zero with increasing wavenumber in a sawtooth power law fashion. 
An example spectrum is shown in Fig. 9. The envelope of the high wavenumber sawtooth 
approximately obeys A;^^. 

If behavior where an exponential high wavenumber tail or a power law with 7 > 2 is 
to be allowed, this type of initial condition should not be used. In fact, Kerr (1993) found 
that 7 > 3 and associated 7 — > 3 as the possible singular time was approached with the test 
that Qpr ~ l/(T-t) 2. 

Increasing the resolution used with the unfiltered initial conditions did not alter either 
the low wavenumber components of the spectra or the magnitude of the regions of negative 
vorticity. That is, regions of negative vorticity in roughly the same location and of the 
same magnitude appear. This type of spectral behavior is reminiscent of Gibbs phenomena, 
that is high wavenumber oscillations that occur when a Galerkin method tries to represent 
physical space discontinuities. This suggests that the origin of the sawtooth spectra and 
negative regions of vorticity is the sharp cutoff in the vorticity profile used by Melander 
& Hussain (1989), although possible inconsistencies between the initial perturbation and 
periodic boundary conditions might also have an influence. This is despite the analytic 
formula where all derivatives smoothly went to zero at the edge of the vortex tube. While 
this might work analytically, all that a numerical calculation sees is that the vorticity goes 
from a finite value on one side of a boundary to continuously zero on the other side. 

5.3 Effect when there is no initial filter 

The calculation in this section is meant to show that if the high wavenumber filter is not 
applied to the initial conditions, then the strain saturates, afterwhich the growth of the 
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Figure 8: y vorticity perpendicular in the symmetry plane, t — 0, t — 9 and t — 13 for Kerr 
(1993) calculation. 
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2;-spectra for unfiltered initial condi- 
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Figure 10: Peak strain along the vor- 
ticity (triangle), peak vorticity (cir- 
cle), peak negative vorticity (plus) 
for unfiltered initial conditions vs. 
time. 



peak vorticity is exponential. Qualitatively everything about this calculation is the same 
as the initiahzation used by Kerr (1993). The only significant difference was that the high 
wavenumber hyperviscous filter was not used. 

Only enough resolution to demonstrate these points was used. There are some weak signs 
that at the last time calculated that the strain is beginning to increase and singular behavior 
might be starting. However, given the success of the filtered initial conditions there did not 
seem to be any reason to follow this further. 

The suggestion at the end of the last subsection is that isolated points of negative vor- 
ticity might cause a problem. The best evidence that these negative regions do affect the 
calculation is that while they were initially very small, the order of 4 x 10~^ compared to 
the initial peak vorticity of 1, by the end of the calculation the negative peak was the same 
order of magnitude as the primary vortex. This is shown in Fig. 10. 

Figure 10 also shows the strain, which saturates, leading to exponential growth of peak 
vorticity. It was also found that the position of the peak of a is the same distance from the 
dividing plane (in z) as the peak vorticity, whereas Kerr (1993) found that the strain peak 
is always somewhat further from the dividing plane. 

What are the differences in the structure? The second frame at t = 9 shows that the 
vorticity that develops from initial condition Fig. 8 is not significantly different than the 
initial flattening from Kerr (1993), except it is remaining flat longer. At the slightly later 
time oi t — 13, a head-tail structure is found as in earlier viscous calculations of vortex 
reconnection. The difference with Kerr (1993) is that the highest contour is spread across 
a vortex sheet above the lower dividing plane, across from its mirror image. 

The spreading of this innermost contour is the best indication that similar behavior might 
be happening in other calculations. 

How might these negative regions create this change in structure? Detailed investigation 
of contours in preparing Kerr (1992) showed that the largest negative regions had been 
sucked into the region between the primary vortices, that is along the dividing plane. They 
then grow enormously and pair with their opposite number across the lower dividing plane 
and induce velocity opposite to that from the primary vorticity. This pulls the tail out and 
creates more flattening. 

A comparison with the late times from Pumir & Siggia (1990) is now useful. That 
is for t > 8.1, beyond when ap ~ 0.05||a;||oo . Recall that in Fig 7, the strain was also 
observed to saturate at late times, with exponential growth of the peak vorticity. Based 
on these comparisons, the late times of Pumir & Siggia (1990) should be checked to see if 
there are similar regions of negative vorticity playing a major role, and if so to determine 
whether these regions might then explain the saturation of the strain in their calculations. 
The highest vorticity contour is now spread over the entire tail and the strain contours are 
not as concentrated. 

Figure 12 shows a three-dimensional isosurface plot of the vorticity squared at t = 12 
for the unfiltered initial condition and Fig. 13 an isosurface for the filtered initial condition 
at t — 12. Overall, the three-dimensional structure for the unfiltered initial condition does 
not have two-dimensionalization of the vortices, one of the properties observed by Pumir 
& Siggia (1990), but does show less deformation than the filtered case from Kerr (1993). 
For the filtered initial condition, a crease in Fig. 13 where the isosurface moves away from the 
dividing plane indicates where three-dimensional curvature of the vortex lines is occurring. 
A similar crease does not appear for the unfiltered case. 



Figure 11: ||w||oo and a contours in the symmetry plane for t — 8.3 from Pumir & Siggia 
(1990). 




Figure 12: Three-dimensional isosurface plot Figure 13: Three-dimensional isosurface plot 

of the vorticity squared for (x, y, z) = (4.9 : of the vorticity squared at i = 12 for the fil- 

8.8,0 : 3.9,0 : 3.9) at i = 12 for unfiltered tered initial conditions in Kerr (1993). 
initial conditions. 



Figure 14: Contour plots through the symmetry plane from Hou & Li (2006) at t = 15 and 
17 (left and right) with the same factor of 5 expansion in the vertical z direction. The t — Yl 
frame has been blown up by a factor of 6. Unlike Fig. 1, the t = 15 contours here look more 
like the t = 6 contours in Fig 5. Or t = 12, except there is no upward tilt on the left, even 
some downward tilt, t = 17 is very similar to the calculation with initial noise in Fig. 5 

6 Discussion of new Hou & Li (2006) calculation 

Hou & Li (2006) have recently released results where the initial condition of Kerr (1993) was 
used exactly as originally programmed by using my old Fortran coding. That is: the same 
initial core profile, the same initial perturbation, and the same high-wavenumber filtering. 
As a result, exactly the same initial behavior is expected and seems to be followed. 

The numerical method was similar to that originally used by Kerr & Hussain (1989). 
That is pseudospectral using sine and cosine transforms to apply the symmetries and different 
resolution in the different directions. The ratios of resolutions used was different than in Kerr 
& Hussain (1989), but not significantly different. The highest resolution was 1536 x 1024 x 
3072. 3072 mesh points in the direction of greatest collapse makes the local resolution at 
the dividing plane comparable to adaptive Chebyshev method of Kerr (1993), while in the 
other directions the maximum resolutions were respectively 3 and 4 times more than those 
used by the most resolved case of Kerr (1993). It was acknowledged by Kerr (2005) that 
the resolution in x was inadequate, although a simple post-processing filter did allow new 
analysis in Kerr (2005). 

Another major difference was how the high wavenumbers were truncated. This report 
will conclude that this is the most significant difference. Both the hybrid pseudospec- 
tral/Chebyshev code reported in Kerr (1993) and the totally pseudospectral results from 
around 1991 discussed here used the standard 2/3rds dealiasing, which turns a pseudospec- 
tral calculation into a true spectral calculation. Hou & Li (2006) claim that one can attain 
better convergence by allowing energy in wavenumbers higher than the 2/3rds limit within 
a prescribed envelope. The dangers of doing this are discussed here in section 4. 

The comparison of the results of Hou & Li (2006) and Kerr (1993) can be divided into 
the following three stages. 

• Early times, t < 10, where no results are reported by Hou & Li (2006) and for which 
this report will assume has identical behavior to Kerr (1993). 
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Figure 15: Comparison of enstrophy and cn- 
strophy production between Hou & Li (2006) 
and Kerr (1993) assuming that they have the 
same values at t — 10. The last two frames 
are fip^ not Q as indicated. In Hou & Li 
(2006) the growth of enstrophy is greater for 
12 < t < 17, but is less for t > 17. The best 
signature of this is the last frame of flpr- 



Intermediate times 10 < t < 16.5 where differences are noted here but not by Hou & Li 
(2006). 



• Late times where t > 16.5 for which differences are reported by Hou & Li (2006) and 
where they claim that they can successfully calculate until t — 19. 

Resolution and late times will be considered first so as to review the differences noted by 
Hou & Li (2006). Then the intermediate stage will be discussed to show that differences 
appear much earlier. This will be used to understand the late time differences. 

6.1 Resolution and late times 

Besides the calculation discussed in Sec. 5.3 that used a fully pseudospectral calculation, 
among the tests done prior to publication of the results in Kerr (1993) was a fully pseu- 
dospectral calculation using nearly the same filtered, initial condition as Kerr (1993). The 
same initial peak vorticity and perturbation were used. The only difference was a less severe 
initial filter. The maximum resolution used was 256 x 64 x 512 and the last reliable time 
was estimated to be i = 15 when ||w||oo — 5.7. This calculation is mentioned here only to 
note that tests were done to determine what level of resolution would be needed to reach 
a given level of singular behavior with a totally pseudospectral code. My recollection that 
a remesh onto the 256 x 64 x 512 mesh from a 128 x 64 x 256 was required at t = 13.5 
when Halloo =4. That is for every desired doubling in ||u?||oo , the resolution needed to be 
quadrupoled in at least the x and z directions. Based upon this estimate and assuming the 
singular growth reported by Kerr (1993), with 3072 mesh points in z, the maximum ||c«;||oo 
that could be obtained would be 14 at t ~ 17.3. 

Based upon this, the claim by Hou & Li (2006) that they can calculate until t = 19 
would only be possible if there isn't singular growth in ||a;||oo and the position of ||u;||oo does 
not collapse as much as in Kerr (1993). As evidence that their calculation is resolved until 
that time, Hou & Li (2006) state that they still have 8 mesh points between the position 
of II^A^Iloo and the dividing plane at that time. A useful comparison between Kerr (1993) 
and Hou & Li (2006) is Fig. 9 in Hou & Li (2006), where there is a clear break in ||w||oo 
starting at t = 16.5. from the ||w||oo ~ c/ (18.7 — t) result of Kerr (1993). 

Prom this one might conclude that the calculation of Kerr (1993) should not have been 
run beyond about t — 16.5. However, the structural differences arise far before t — 16.5. 

6.2 Structural differences at t = 15 

Based on the comparisons of resolutions in Kerr (1993) and the pseudospectral test just 
mentioned, the evidence is that Kerr (1993) is resolved at least until t = 16.5. Fig. 1, 
take from Kerr (1993) at t = 15 < 16.5, shows that the innermost vorticity contour is in 
the corner. In Fig. 14, taken from Fig. 15 of Hou & Li (2006), the innermost contour is 
flattened. 

What the contours from Hou & Li (2006) most resemble are the later times from the 
unfiltcrcd calculation from Kerr (1992) in Fig. 8 and t > 8.2 from Pumir & Siggia (1990) 
shown in Fig. 11. This suggests that noise at the small scales that becomes amplified is the 
source of the differences. Since they do filter their initial condition, where does the necessary 
noise come from? The most likely source would be retaining wavenumbers beyond the 2/3rds 
rule cutoff. 



The simplest test of this proposal would be to check for small negative regions of vorticity, 
locate them, and plot how they grow in time. Our plan is to run independent calculations 
to check this possibility. 

Another sign of a difference for t < 17 is obtained by comparing the behavior of enstrophy 
and its production plotted in Fig 5.3. Upon my suggestion, Hou & Li (2006) added the time 
dependence of enstrophy Q and enstrophy production fip^ to the version they submitted. 
The units arc not those used by Kerr (1993), so in Fig 5.3 the assumption is made that 
the Kerr (1993) and Hou & Li (2006) calculations are identical for t < 10 and therefore 
by scaling their values at t = 10, a comparison can be made. The discussion in Hou & Li 
(2006) suggests that the values of ilpr shown were calculated from {d/dt)Vt and not from 
J dVcUiCijCUj directly. 

Both Q and Qpr grow faster than in Kerr (1993) for 12 < t < 16 and then more slowly. 
This is seen best in the third figure for Qpr. This trend would be consistent with there being 
more enstrophy growth in the long tail initially, but less concentration in the corner to build 
upon at later times. In the table below the numbers used to make these plots are given. The 
enstrophy results from Kerr (1993) correct a misprint. The formula for enstrophy in Kerr 
(1993) should have been: ft = -.0105 log(18.9 - t) + .05, giving Q{t = 0) = .019. 
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7 Moving forward 

In March 2003, A. Bhattacharee, U. Frisch, R.M. Kerr, N. Zabusky and others met at the 
Institute for Advanced Studies in Princeton to consider what direction a new computational 
effort should follow. This meeting was instigated by the untimely death of our friend and 
coUegue Rich Pelz. 

The primary conclusion that was reached was that there was a serious problem where 
each numerical method seemed to give a different answer. It was concluded that this was 
primarily because each team was pushing their calculations a little too far. In the early 
1990s there was really no choice. To get any trend one had to overextend the calculations. 
In the first decade of the 21st century there is the luxury of far greater computing power 
and significantly improved numerical methods. Therefore it was decided to propose an 
international collaborative effort based on: 

• A factor of at least 4 increase in resolution in each direction 

• Availability of modern adaptive mesh finite difference and spectral element codes ca- 
pable of providing significant improvements in local resolution. 

• Separate groups should conduct simultaneous simulations with the same initial con- 
ditions. The only results that should be reported are until a time when the two 
calculations differ only by an agreed upon amount. 

Only in the UK was a proposal made. A separate proposal with no communication with 
the people above, and therefore not informed of the rules, was funded and went to T. Hou 



at Cal Tech, R. Caflisch at UCLA and M. Siegel at the New Jersey Institute of Technology. 
The UK proposal between Warwick and Imperial was never funded due to the poor rules for 
reviewing interdisciplinary proposals in the UK. Funding has finally been obtained from the 
Leverhulme Foundation, but only for one post-doc at Warwick. 

The original intention was that the two first adaptive codes to be used would be the 
spectral element code of Barkley at Warwick and some variant of the finite-difference adaptive 
code of Graucr in Bochum or from within Bhattacharjeee's group at New Hampshire. Later, 
the improved spectral element code of Sherwin from Imperial would be used. 

Most likely at this time we will go directly to Sherwin's spectral element code and work 
with Pumir on reviving his adaptive-mesh finite difference code. 

However, the first priority now is to address the issues raised by Hou & Li (2006). There- 
fore, their calculations will be repeated and compared in detail with a straight pseudospectral 
calculation using a strict 2/3rds rule and stopped early enough so that there are no questions 
about small-scale resolution. This should be near t = 17, although the comparison above 
suggests that a difference should show up as early as i = 15. 
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